/*--------------实验三---------------*/
# include <stdio.h>
# include <math.h>

# define PI 3.14
# define G 6.67e-5
# define Q 2.67

/*-----函数声明-----*/
double Pan_d(double b, double a);
double R(double x, double y, double z, double x_k, double y_k, double z_k);
double San_G(double x, double y, double z, double x_k, double y_k, double z_k);
double San_X(double x, double y, double z, double x_k, double y_k, double z_k, double M_x, double M_y, double M_z);
double San_Y(double x, double y, double z, double x_k, double y_k, double z_k, double M_x, double M_y, double M_z);
double San_Z(double x, double y, double z, double x_k, double y_k, double z_k, double M_x, double M_y, double M_z);

void main()
{
	int i, j, k, l, p;//循环变量部分 
	double x_o=1000.0, y_o=1000.0, z_o=1000.0, a=200.0, b=200.0, c=200.0, M_x=200.0, M_y=200.0, M_z=200.0;
	double x_k[101], y_k[101], z_k=0, sum_g=0, sum_x=0, sum_y=0, sum_z=0, sAn_g[101][101], sAn_x[101][101], sAn_y[101][101], sAn_z[101][101];
	/*--------方便带入稍后解释-------*/
	double X_ju[2]={x_o-(a/2), x_o+(a/2)}, Y_ju[2]={y_o-(b/2), y_o+(b/2)}, Z_ju[2]={z_o-(c/2), z_o+(c/2)};
	FILE *fp; //文件指针
	fp = fopen("BFHBFHBFH.txt", "w+");
	
	/*------算法部分------*/ 
	for(i=0; i<101; i++) {
		/*----101*101点赋值------*/
		x_k[i] = i*(20.0);
		y_k[i] = i*(20.0); 
		for(j=0; j<101; j++) {
			/*-----8次循环固定算法（积分部分）----*/
			for(k=0; k<2; k++) {
				for(l=0; l<2; l++) {
					for(p=0; p<2; p++) {
						sum_g += (pow(-1, (p+k+l)) * San_G(X_ju[k], Y_ju[l], Z_ju[p], x_k[i], y_k[j], z_k));
						sum_x += (pow(-1, (p+k+l)) * San_X(X_ju[k], Y_ju[l], Z_ju[p], x_k[i], y_k[j], z_k, M_x, M_y, M_z));
						sum_y += (pow(-1, (p+k+l)) * San_Y(X_ju[k], Y_ju[l], Z_ju[p], x_k[i], y_k[j], z_k, M_x, M_y, M_z));
						sum_z += (pow(-1, (p+k+l)) * San_Z(X_ju[k], Y_ju[l], Z_ju[p], x_k[i], y_k[j], z_k, M_x, M_y, M_z));
					}	
				}
			}
			
			sAn_g[i][j]	= sum_g;
			sAn_x[i][j] = sum_x;
			sAn_y[i][j] = sum_y;
			sAn_z[i][j] = sum_z;
			/*------写入文件------*/
			printf("[%lf][%lf]\n", x_k[i], y_k[j]); 
			fprintf(fp, "%lf %lf %lf %lf\n", sAn_g[i][j], sAn_x[i][j], sAn_y[i][j], sAn_z[i][j]);
		}		
	}
	/*------关闭文件------*/
	fclose(fp);
}

double Pan_d(double b, double a)
{
	if(fabs(a)>1.0e-15) {
		if(b/a > 0.0) {
			if(a>0.0 && b>0.0) {
				return atan2(b,a);
			}
			else {
				return (-PI + atan2(b,a));
			}
		}
		else {
			return (PI + atan2(b,a));
		}
	}
	else {
		if(b<0.0) {
			return -PI/2;
		}
		else {
			return PI/2;
		}
	}
}

double R(double x, double y, double z, double x_k, double y_k, double z_k)
{
	return sqrt((pow(x-x_k, 2) + pow(y-y_k, 2) + pow(z-z_k, 2)));
}

double San_G(double x, double y, double z, double x_k, double y_k, double z_k)
{	
	return (-1*G*Q*((x-x_k)*log(y-y_k+R(x,y,z,x_k,y_k,z_k)) + (y-y_k)*log(x-x_k+R(x,y,z,x_k,y_k,z_k)) + ((z-z_k)*Pan_d((x-x_k)*(y-y_k),(z-z_k)*R(x,y,z,x_k,y_k,z_k) ) )));
}

double San_X(double x, double y, double z, double x_k, double y_k, double z_k, double M_x, double M_y, double M_z)
{
	return ((1/(4*PI))*((-1)*M_x*Pan_d((y-y_k)*(z-z_k),R(x,y,z,x_k,y_k,z_k)*(x-x_k)) + M_y*log(R(x,y,z,x_k,y_k,z_k)+(z-z_k))+M_z*log(R(x,y,z,x_k,y_k,z_k)+(y-y_k))));
}

double San_Y(double x, double y, double z, double x_k, double y_k, double z_k, double M_x, double M_y, double M_z)
{
	return ((1/(4*PI))*((-1)*M_y*Pan_d((x-x_k)*(z-z_k),R(x,y,z,x_k,y_k,z_k)*(y-y_k)) + M_x*log(R(x,y,z,x_k,y_k,z_k)+(z-z_k))+M_z*log(R(x,y,z,x_k,y_k,z_k)+(x-x_k))));
}

double San_Z(double x, double y, double z, double x_k, double y_k, double z_k, double M_x, double M_y, double M_z)
{
	return ((1/(4*PI))*((-1)*M_z*Pan_d((x-x_k),(z-z_k))*((y-y_k)/R(x,y,z,x_k,y_k,z_k)) + M_x*log(R(x,y,z,x_k,y_k,z_k)+(y-y_k))+(-1)*M_y*log(R(x,y,z,x_k,y_k,z_k)+(x-x_k))));
}

